﻿using System;
using System.Text;
using System.Drawing;
using System.Buffers;
using System.Collections;
using System.Collections.Generic;
using System.Runtime.InteropServices;

public static partial class NativeAOT
{
    [UnmanagedCallersOnly(EntryPoint = "pqshoot")]
    public static unsafe double pqshoot(int n, double a, double b, double eps, IntPtr y_ptr, IntPtr f_abc_ptr)
    {
        double* y = (double*)y_ptr.ToPointer();
        f_abc = Marshal.GetDelegateForFunctionPointer<delegatefunc_abc>(f_abc_ptr);

        return pqshoot(n, a, b, eps, y);
    }

    /// <summary>
    /// 求解二阶边值问题连分式法
    /// f计算二阶常微分方程右端函数f(t,y,z)
    /// </summary>
    /// <param name="n">积分区间等分数</param>
    /// <param name="a">积分区间左端点</param>
    /// <param name="b">积分区间右端点。要求b>a。</param>
    /// <param name="eps">控制精度要求</param>
    /// <param name="y">y[n+1]y[0]存放左端点边界值y(a)，y[n]存放右端点边界值y(b)。返回n+1个等距离散点上的数值解。</param>
    /// <returns>函数返回y在左端点处的一阶导数值</returns>
    public static unsafe double pqshoot(int n, double a, double b, double eps, double* y)
    {
        int i, j, il, flag;
        double y0, z0, t, h;
        double* bb = stackalloc double[10];
        double* zz = stackalloc double[10];
        double* yn = stackalloc double[10];

        h = (b - a) / n;
        il = 0; z0 = 0.0;
        flag = 0;
        while ((il < 20) && (flag == 0))
        {
            il = il + 1;
            j = 0;
            zz[0] = z0;
            t = a;
            y0 = y[0];
            // 计算yn[0]
            for (i = 1; i <= n; i++)
            {
                pqeuler2(t, h, &y0, &z0, eps);
                t = t + h;
            }
            yn[0] = y0;
            // 计算bb[0]
            bb[0] = zz[0];
            j = 1;
            zz[1] = zz[0] + 0.1;
            z0 = zz[1];
            t = a;
            y0 = y[0];
            // 计算yn[1]
            for (i = 1; i <= n; i++)
            {
                pqeuler2(t, h, &y0, &z0, eps);
                t = t + h;
            }
            yn[1] = y0;
            while (j <= 7)
            {
                //计算bb[j]
                funpqj(yn, zz, bb, j);
                // 计算zz[j+1]
                zz[j + 1] = funpqv(yn, bb, j, y[n]);
                z0 = zz[j + 1];
                t = a;
                y0 = y[0];
                // 计算yn[j+1]
                for (i = 1; i <= n; i++)
                {
                    pqeuler2(t, h, &y0, &z0, eps);
                    if (i < n) y[i] = y0;
                    t = t + h;
                }
                yn[j + 1] = y0;
                z0 = zz[j + 1];
                if (Math.Abs(yn[j + 1] - y[n]) >= eps) j = j + 1;
                else j = 10;
            }
            if (j == 10)
            {
                flag = 1;
            }
        }
        return (z0);
    }

    /*
    // 求解二阶边值问题连分式法例1
        int main()
        {
            int k;
            double dy0, y[11], f(double,double,double);
            y[0] = 0.0;  y[10] = 1.0;
            dy0 = pqshoot(10, 0.0, 1.0, 0.0000001, y, f);
            cout <<"初始斜率 = " <<dy0 <<endl;
            for (k=0; k<11; k++)
                cout <<"x = " <<0.1*k <<"     y = " << y[k] <<endl;
            return 0;
        }
    // 计算二阶微分方程右端函数值
        double  f(double t, double y, double z)
        {
            double d;
            d = t + y;
            return(d);
        }
    */
    /*
    // 求解二阶边值问题连分式法例2
        int main()
        {
            int k;
            double dy0, y[11], f(double,double,double);
            y[0] = 1.0;  y[10] = 2.0;
            dy0 = pqshoot(10, 0.0, 1.0, 0.0000001, y, f);
            cout <<"初始斜率 = " <<dy0 <<endl;
            for (k=0; k<11; k++)
                cout <<"x = " <<0.1*k <<"     y = " << y[k] <<endl;
            return 0;
        }
    // 计算二阶微分方程右端函数值
        double  f(double t, double y, double z)
        {
            double d;
            d = (6*t-3.0+t*z+3*y)/(1.0+t*t);
            return(d);
        }
    */
}

